# 进行学习连续型变量学习 2023年3月3日14:31:37 
#### 0 准备好环境####
library(tidyverse)
library(gtsummary)
library(haven)
library(tidyr)
library(survey)

#### 1 常规线性模型 ####
##### 1.1 生成数据 #####
n = 30
#rnorm 正态分布 参数 数量 均值 标准差
# 生成体重数据
weight.male <- rnorm(n,75,15)
weight.female <- rnorm(n,55,10)
weight <- c(weight.male,weight.female)
# sample() 的含义是随机把数组进行打乱排序  
# 生成性别变量
sex.male <-factor(sample(c('male'),n,T))
sex.famale <-factor(sample(c('famale'),n,T))
sex <- c(sex.male,sex.famale)

# 生成体重因素
height <- weight + rnorm(n *2 ,0,5)
##### 1.2 生成演示用的数据信息 #####

# 构建简单的线性回归模型
data <- data.frame(height = height, weight = weight, sex = sex)
# 使用glm 进行简单回归模型
lm <- glm(weight ~ sex+height,data = data) # y ~ x1+x2 
lm2 <- glm(height ~ sex+weight,data = data) # y ~ x1+x2 


# 大体看下数据 模型详细结果
summary(lm) 
# 列出模型参数 系数
coefficients(lm)
# 提供模型的置信区间 
confint(lm)
# 列出模型aic 信息统计量??? 什么意思 这个值越低越好点
AIC(lm)
AIC(lm2)
# 用拟合模型对y 值的预测? 什么意思
predict(lm)

##### 1.2 使用 gtsummary 优化输出 #####
# 生成了论文类似的数据
lm_tbl <- tbl_regression(lm)

#### 2. 加权线性模型 ####
##### 2.1 提取数据模块 #####
setwd('G:/BaiduNetdiskDownload/NHANES')
demo.i <- read_xpt("2015-2016/Demographics/demo_i.xpt")#要提前设置好数据存储的路径
colnames(demo.i)
analyze.sample.data <- demo.i[,c('SEQN',
                                 'RIDRETH3', "RIDAGEYR", "RIAGENDR", "DMDEDUC2", "INDFMPIR", 
                                 "WTINT2YR", "SDMVPSU", "SDMVSTRA")]

# 一键清除na的行
analyze.sample.data.drop.na <- drop_na(analyze.sample.data) 
dim(analyze.sample.data.drop.na) 
# AD <- analyze.sample.data[which(analyze.sample.data[,"SEQN"]==83738),]
# print(AD)
# AD2 <- analyze.sample.data.drop.na[which(analyze.sample.data.drop.na[,"SEQN"]==83738),]
# print(AD2)

# 设计复杂抽样包 默认执行 因为是两年 权重就取  WTINT2YR
NHANES_design <- svydesign(
  data = analyze.sample.data.drop.na, 
  ids = ~SDMVPSU, 
  strata = ~SDMVSTRA, 
  nest = TRUE, 
  weights = ~WTINT2YR
)
##### 2.4 线性回顾模型讲解 #####
# 年龄、教育和贫困指数之间的相关性
lm <- svyglm(INDFMPIR ~ RIDAGEYR + factor(DMDEDUC2), design = NHANES_design)


# 通用函数 
summary(lm) 
# 列出模型参数 系数
coefficients(lm)
# 提供模型的置信区间 
confint(lm)
# 列出模型aic 信息统计量??? 什么意思 这个值越低越好点
AIC(lm)
# 用拟合模型对y 值的预测? 什么意思
predict(lm)

##### 1.2 使用 gtsummary 优化输出 #####
# 生成了论文类似的数据
lm_tbl <- tbl_regression(lm)

tbl_regression(lm)%>% 
  add_global_p() %>%
  
  add_glance_table(
    # include = c(nobs, AIC, BIC))
    include = c( nobs, AIC, BIC))


